X=imread('wed2.jpg');

%show original picture 600x600
Xoriginal=X(1:600,1:600,:);
f1=figure;
set(f, 'name', 'Original 600x600','numbertitle','off')
imagesc(Xoriginal);
axis square;

r=X(:,:,1);g=X(:,:,2);b=X(:,:,3); %seperate colors

for i=1:75
    for j=1:75
        xr=r((i-1)*8+1:(i-1)*8+8,(j-1)*8+1:(j-1)*8+8); %extract 8x8 pixel block
        xg=g((i-1)*8+1:(i-1)*8+8,(j-1)*8+1:(j-1)*8+8); %extract 8x8 pixel block
        xb=b((i-1)*8+1:(i-1)*8+8,(j-1)*8+1:(j-1)*8+8); %extract 8x8 pixel block

        %Start of Compression Process
        xrd=double(xr); %convert uint8 to double
        xgd=double(xg); %convert uint8 to double
        xbd=double(xb); %convert uint8 to double

        xrc=xrd-128; %center matrix over zero
        xgc=xgd-128; %center matrix over zero
        xbc=xbd-128; %center matrix over zero

        Yr=dct(dct(xrc')'); %Apply the 2D-DCT, then Y is the transform matrix
        Yg=dct(dct(xgc')'); %Apply the 2D-DCT, then Y is the transform matrix
        Yb=dct(dct(xbc')'); %Apply the 2D-DCT, then Y is the transform matrix

        p=1;    %p is the loss parameter
        Q=p*8./hilb(8); %use hilbert matrix to define the linear quantization matrix
        Yrq=round(Yr./Q); %replace Y by the compressed matrix (Quantize)
        Ygq=round(Yg./Q); %replace Y by the compressed matrix (Quantize)
        Ybq=round(Yb./Q); %replace Y by the compressed matrix (Quantize)

        %Start of Decompression Process
        Yrdq=Yrq.*Q; %dequantization
        Ygdq=Ygq.*Q; %dequantization
        Ybdq=Ybq.*Q; %dequantization

        Xrdq=idct2(Yrdq); %inverse DCT transform
        Xgdq=idct2(Ygdq); %inverse DCT transform
        Xbdq=idct2(Ybdq); %inverse DCT transform

        Xre=Xrdq+128; %un-center
        Xge=Xgdq+128; %un-center
        Xbe=Xbdq+128; %un-center

        Xout((i-1)*8+1:(i-1)*8+8,(j-1)*8+1:(j-1)*8+8, 1)=uint8(Xre); %convert back to uint8
        Xout((i-1)*8+1:(i-1)*8+8,(j-1)*8+1:(j-1)*8+8, 2)=uint8(Xge); %convert back to uint8
        Xout((i-1)*8+1:(i-1)*8+8,(j-1)*8+1:(j-1)*8+8, 3)=uint8(Xbe); %convert back to uint8

    end
end

f2=figure;
set(f2,'name', 'After, p=1','numbertitle','off');
imagesc(Xout); %display the image after compression
axis square;

f3=figure;
set(f3,'name', 'Comparison with p=1','numbertitle','off');
h = imshowpair(Xoriginal,Xout,'diff');
axis square;